#ifndef EIGEN_USE_MKL_ALL
#define EIGEN_USE_MKL_ALL
#endif

#include <Eigen/Dense>
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <map>
#include <set>
#include <numeric>
#include <ctime>

using namespace Eigen;
using namespace std;

const int times = 10000;

namespace cart_tree
{
	template <typename T>
	struct ColComparator
	{
		const MatrixBase<T>& matrix;
		const int row;

		ColComparator(MatrixBase<T>& m, int row)
			: matrix(m.derived()), row(row)
		{
		}

		bool operator()(int i, int j) const
		{
			return matrix(row, i) < matrix(row, j);
		}
	};

	const double MIN_DOUBLE = 1e-9;

	struct MyCompDouble
	{
		bool operator()(const double& left, const double& right) const
		{
			double x = left - right;
			if (x < 0 && abs(x) > MIN_DOUBLE)
				return true;
			return false;
		}
	};


	//need verify
	void sortByCol(MatrixXd& m, int n)
	{
		MatrixXd tmp = m;
		VectorXi indice = VectorXi::LinSpaced(m.cols(), 0, m.cols() - 1);
		std::sort(indice.data(), indice.data() + indice.size(),
			ColComparator<MatrixXd>(m, n));
		for (int i = 0; i < indice.size(); i++)
			m.col(i) = tmp.col(indice(i));
	}

	//need verify
	void premuteByCol(MatrixXd& m)
	{
		MatrixXd tmp = m;
		VectorXi indice = VectorXi::LinSpaced(m.cols(), 0, m.cols() - 1);
		std::random_shuffle(indice.data(), indice.data() + indice.size());
		for (int i = 0; i < indice.size(); i++)
			m.col(i) = tmp.col(indice(i));
	}

	MatrixXd* readFromFile(const string& fileName)
	{
		ifstream ifd(fileName, ifstream::in);
		string line;
		stringstream ss;
		int rows, cols;
		if (!ifd.is_open())
		{
			cerr << "Open Error: " << fileName;
			exit(1);
		}
		getline(ifd, line);
		ss << line;
		ss >> rows >> cols;
		ss.clear();
		MatrixXd* m = new MatrixXd(rows, cols);
		int c = 0;
		while (std::getline(ifd, line))
		{
			ss << line;
			for (int i = 0; i < cols; ++i)
			{
				ss >> (*m)(c, i);
			}
			ss.clear();
			++c;
		}
		ifd.close();
		m->transposeInPlace();
		return m;
	}

	//last row is the label
	double calcGini(const MatrixXd& m)
	{
		ArrayXd lbs = m.row(m.rows() - 1);
		ArrayXd prop(lbs.maxCoeff());
		for (int i = 1; i <= lbs.maxCoeff(); i++)
			prop(i - 1) = (lbs == i).count();
		prop = prop / (prop.sum());
		return 1 - prop.matrix().squaredNorm();
	}

	void binSplit(const MatrixXd& m, int iFeature, double vFeature, MatrixXd& left, MatrixXd& right)
	{
		int lc = (m.row(iFeature).array() < vFeature).count();
		int rc = m.cols() - lc;
		if (lc)
			left.resize(m.rows(), lc);
		if (rc)
			right.resize(m.rows(), rc);
		int li = 0, ri = 0;
		for (int i = 0; i < m.cols(); i++)
		{
			if (m(iFeature, i) < vFeature)
				left.col(li++) = m.col(i);
			else
				right.col(ri++) = m.col(i);
		}
	}

	//for getting intervals of set v
	template <typename T, typename U>
	VectorXd getIntervals(const set<T, U>& v)
	{
		VectorXd r(v.size() - 1);
		vector<double> tmp;
		for (auto p = v.cbegin(); p != v.cend(); ++p)
			tmp.push_back(*p);
		for (int i = 0; i < tmp.size() - 1; ++i)
			r(i) = (tmp[i] + tmp[i + 1]) / 2;

		return r;
	}

	//for getting different values for Matrix m at i-th Feature
	//T for element type,U for Matrix Type, V for compare
	template <typename T, typename U, typename V>
	void getDiffValues(const MatrixBase<U>& m, int iFeature, set<T, V>& s)
	{
		s.clear();
		for (int i = 0; i < m.cols(); ++i)
			s.insert(m(iFeature, i));
	}


	void chooseBestSplit(const MatrixXd& m, int& iFeature, double& vFeature, const vector<int> &flags)
	{
		set<double, MyCompDouble> fvs;
		double gini = calcGini(m), giniTmp;
		double maxInfo = (numeric_limits<double>::min)();
		iFeature = -2;//if every feature value is equal, return -2 to let the parent know
		vFeature = 0;
		//i-th feature
		for (int i = 0; i < m.rows() - 1; ++i)
		{
			if (flags[i])
				continue;
			getDiffValues<double, MatrixXd, MyCompDouble>(m, i, fvs);//if all vFeature for i-th Feature is same, the size of fvs will be one, and no intervals
			VectorXd intervals = getIntervals(fvs);
			for (int j = 0; j < intervals.size(); ++j)
			{
				MatrixXd left, right;
				binSplit(m, i, intervals[j], left, right);
				giniTmp = (calcGini(left) * left.size() + calcGini(right) * right.size()) / m.size();
				giniTmp = gini - giniTmp;
				if (giniTmp > maxInfo)
				{
					maxInfo = giniTmp;
					iFeature = i;
					vFeature = intervals[j];
				}
			}
		}
	}

	struct Cart_Tree
	{
		int iFeature;
		double vFeature;
		int label;
		MatrixXd leafs;
		Cart_Tree *left, *right;
	};

	//so diao
	template<typename T>
	bool pairCompare(const T &l, const T &r)
	{
		return l.second < r.second;
	}
	template<class T>
	typename T::const_iterator map_max_element(const T &A)
	{
		typedef typename T::value_type pair_type;
		return max_element(A.begin(), A.end(), pairCompare<pair_type>);
	}


	//create tree, train
	Cart_Tree* createTree(const MatrixXd& m, vector<int> flags)
	{
		if (m.size() == 0)
			return nullptr;
		map<int, int> mm;
		for (auto i = 0; i < m.cols(); ++i)
			mm[static_cast<int>(m(m.rows() - 1, i))]++;
		//if all labels in m are same return tree
		if (mm.size() == 1)
		{
			Cart_Tree* tree = new Cart_Tree;
			map<int, int>::const_iterator x = mm.cbegin();
			tree->iFeature = -1;
			tree->label = x->first;
			tree->left = nullptr;
			tree->right = nullptr;
			tree->leafs = m;
			return tree;
		}
		//if all feature were processed create the leaf
		if (accumulate(flags.cbegin(), flags.cend(), 0) == flags.size())
		{
			Cart_Tree *t = new Cart_Tree;
			map<int, int>::const_iterator x = map_max_element(mm);
			t->iFeature = -1;
			t->label = x->first;
			t->left = nullptr;
			t->right = nullptr;
			t->leafs = m;
			return t;
		}
		int iFeature;
		double vFeature;
		chooseBestSplit(m, iFeature, vFeature, flags);
		//if all vFeature is same create the leaf
		if (iFeature == -2)
		{
			Cart_Tree *t = new Cart_Tree;
			map<int, int>::const_iterator x = map_max_element(mm);
			t->iFeature = -1;
			t->label = x->first;
			t->left = nullptr;
			t->right = nullptr;
			t->leafs = m;
			return t;
		}
		flags[iFeature] = 1;
		MatrixXd left, right;
		binSplit(m, iFeature, vFeature, left, right);
		Cart_Tree* tree = new Cart_Tree;
		tree->iFeature = iFeature;
		tree->vFeature = vFeature;
		tree->label = 0;
		tree->left = createTree(left, flags);
		tree->right = createTree(right, flags);
		return tree;
	}


	//test 
	int getLabel(const Cart_Tree &t, const VectorXd &data)
	{
		if (t.iFeature == -1)
			return t.label;
		if (data(t.iFeature) < t.vFeature)
			return getLabel(*t.left, data);
		return getLabel(*t.right, data);
	}
	double validate(const Cart_Tree &t, const MatrixXd &test)
	{
		int count = 0;
		for (auto i = 0; i < test.cols(); ++i)
		{
			bool tmp = (getLabel(t, test.col(i)) == (int)test(test.rows() - 1, i));
			count += tmp ? 1 : 0;
		}
		return (double)count / test.cols();
	}

	void printTree(const Cart_Tree *t, int count)
	{
		if (t)
		{
			for (int i = 0; i < count; ++i)
				cout << "\t";
			cout << t->iFeature << "," << t->label << endl;
			printTree(t->left, count + 1);
			printTree(t->right, count + 1);
		}
	}
	void eraseTree(Cart_Tree *t)
	{
		if (t->left)
			eraseTree(t->left);
		if (t->right)
			eraseTree(t->right);
		delete(t);
	}

	void getLeaves(const Cart_Tree *t, MatrixXd &m)
	{
		if (t&&t->iFeature == -1)
		{
			MatrixXd tmp = m;
			m.resize(tmp.rows() ? tmp.rows() : t->leafs.rows(), tmp.cols() + t->leafs.cols());
			if (tmp.rows())
				m << tmp, t->leafs;
			else
				m << t->leafs;
			return;
		}
		getLeaves(t->left, m);
		getLeaves(t->right, m);
	}

	//post Prune
	Cart_Tree *postPruning(Cart_Tree *t, const Cart_Tree& root, const MatrixXd &test, const double before)
	{
		if (!t->left && !t->right)
			return t;
		Cart_Tree *left = nullptr, *right = nullptr;
		if (t->left)
			left = postPruning(t->left, root, test, before);
		if (t->right)
			right = postPruning(t->right, root, test, before);
		map<int, int> lbs;
		MatrixXd leaves;
		getLeaves(t, leaves);
		for (int i = 0; i < leaves.cols(); ++i)
			lbs[leaves(leaves.rows() - 1, i)]++;
		auto it = map_max_element(lbs);
		int tmpIfeature = t->iFeature;
		t->iFeature = -1;
		t->label = it->first;
		double tmpRec = validate(root, test);
		if (tmpRec < before)
			t->iFeature = tmpIfeature;
		else
		{
			t->left = nullptr;
			t->right = nullptr;
			eraseTree(left);
			eraseTree(right);
			t->leafs = leaves;
		}
		return t;
	}
}


int main()
{
	MatrixXd* d;
	std::srand(unsigned(std::time(0)));
	d = cart_tree::readFromFile("iris.data");
	double rec = 0;
	for (int i = 0; i < times; i++)
	{
		cart_tree::premuteByCol(*d);
		MatrixXd train, verify, test;
		train = d->block(0, 0, d->rows(), 50);
		verify = d->block(0, 50, d->rows(), 50);
		test = d->block(0, 100, d->rows(), 50);

		vector<int> flags;
		for (int i = 0; i < d->rows() - 1; ++i)
			flags.push_back(0);
		cart_tree::Cart_Tree* t = cart_tree::createTree(train, flags);
		postPruning(t, *t, verify, rec);
		rec += cart_tree::validate(*t, test);
		cart_tree::eraseTree(t);
	}
	cout << rec / times * 100 << endl;
	return 0;
}
